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The problem of ballistically controlled annihilation is revisited for general initial velocity distri- 
butions and arbitrary dimension. An analytical derivation of the hierarchy equations obeyed by the 
reduced distributions is given, and a scaling analysis of the corresponding spatially homogeneous 
system is performed. This approach points to the relevance of the non-linear Boltzmann equation for 
dimensions larger than one and provides expressions for the exponents describing the decay of the 
particle density n(t) oc t~^ and the root mean-square velocity v tx t~"' in term of a parameter related 
to the dissipation of kinetic energy. The Boltzmann equation is then solved perturbatively within a 
systematic expansion in Sonine polynomials. Analytical expressions for the exponents £ and 7 are 
obtained in arbitrary dimension as a function of the parameter p characterizing the small velocity 
behavior of the initial velocity distribution. Moreover, the leading non-Gaussian corrections to the 
scaled velocity distribution are computed. These expressions for the scaling exponents are in good 
agreement with the values reported in the literature for continuous velocity distributions in d = 1. 
For the two dimensional case, we implement Monte-Carlo and molecular dynamics simulations that 
turn out to be in excellent agreement with the analytical predictions. 



I. INTRODUCTION 



Ballistically controlled reactions provide simple examples of non-equilibrium systems with complex kinetics. They 
consist of an assembly of particles with a given velocity distribution, moving freely between collisions in a d-dimensional 
space. In the simplest version of these models we consider here, when two particles meet, they instantaneously 
annihilate each other and disappear from the system. Despite its apparent simplicity, this problem is highly non 
trivial and has attracted substantial interest during the past years (l], [| [| ^, ||, [| ^, ||, ^, 10 . The one-dimensional 
case where the particles can only have two velocities ±c has been studied in a pioneering work by Elskcns and Frisch 
. In particular, they proved that for a symmetric initial velocity distribution, the particle density is decreasing, in 
the long time limit, as n (t) oc r« oc t- 1 / 2 . The case of general distributions in dimension d — 1 was discussed by 
Piasecki [|] , who reduced exactly the annihilation dynamics to a single closed equation for the two-particle conditional 
probability. Moreover, it was shown in the particular bimodal situation of a discrete velocity distribution (±c) that 
in one dimension, important dynamical correlations are developing during the time evolution, invalidating mean-held 
or Boltzmann-like approaches. This exact approach was applied to the case of a three velocity distribution by Droz 
et al H , with the result that the dynamical exponents were strongly depending upon the details of the initial velocity 
distribution. 

No analytical solutions could be found for continuous velocity distributions. In one dimension, Ben-Nairn et al.[|], ^ 
have shown that the exponent £ could depend on the behavior near the origin of the initial velocity distribution. This 
problem has been revisited by Rey et al Q . Based on the exact theoretical approach |^|, [j| , a dynamical scaling theory 
was derived and its validity supported by numerical simulations for several velocity distributions. This leads to the 
conjecture that all continuous velocity distributions f(v) that are symmetric, regular, and such that /(0) 7^ are 
attracted, in the long-time regime, towards the same distribution and thus belong to the same universality class. This 



conjecture was reinforced by numerical simulations in two dimensions 10 



The case of a continuous velocity distribution has also been approached recently by Krapivsky and Sire |9|. Starting 
from a Boltzmann equation, they investigated the decay of the particle density n(t) ~ and the root mean-square 
velocity v oc i~ 7 . They derived upper and lower bounds for the exponents as well as their leading expansion in 1/d, 
valid in high dimension. The main question with such an approach concerns the validity of a Boltzmann equation. 
This is not justified in ID and remains an open problem in higher dimensions. 

The purpose of this paper is to give a first principle answer to this type of question. The article is organized as 
follows. In section 0, an original analytical derivation of the equations governing the dynamics of ballistic annihilation 
is given. The hierarchy equations obeyed by the reduced distributions are obtained. It is shown that in the Grad 
limit, the hierarchy formally reduces to a Boltzmann-like form for d > 1. If the initial reduced distributions factorize, 
the whole hierarchy reduces to one non-linear equation. In the second part of section O, a scaling analysis of the 
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exact spatially homogeneous hierarchy is performed. The exponents £ and 7 are showed to depend only on one 
parameter a related to the dissipation of energy. This scaling analysis turns out to be invalid for the case d = 1 
with discrete velocity distributions, but correct in the continuous case. Strong arguments are given in favor of the 
validity of the Boltzmann approach for the case d > 1 in the long time limit. The Boltzmann equation is then solved 
within a systematic approximation based on an expansion in Sonine polynomials (section III). The first non-Gaussian 
corrections to the scaled velocity distribution are computed and predictions for the exponents £ and 7 are explicitly 
worked out as a function of the dimension d and the parameter \x characterizing the small velocity behaviour of 
the initial velocity distribution: [/(v, t — 0) oc Iv^ for |v| — > 0]. These predictions for £ and 7 are asymptotically 
exact for large dimensions, and reproduce the 1/d correction to the mean-field values. In ID, they are in very good 
agreement with the exponents reported in the literature [^) at the Boltzmann level. In 2D, we implement extensive 
Direct Simulation Monte Carlo methods (DSMC) where the non-linear Boltzmann equation is solved, and Molecular 
Dynamics (MD) simulations where the exact equations of motion are integrated (section fv\) . The agreement between 
the MD and DSMC routes confirms the validity of the Boltzmann approach, and the decay exponents measured are in 
exceptionally good agreement with the Sonine prediction. Conclusions are drawn in section |y|. A preliminary account 
of part of the results presented here has been published elsewhere [fol . 



II. EXACT RESULTS 



A. Derivation of the hierarchy 



Let fl be a region of finite measure in R . We denote by 

/4?(ri, Vl , ... ,r k ,v k ;t) (1) 

the probability density for finding at time t exactly k particles within in the states (r^, Vj) E ft, j — 1,2,..., A;, 
where ya and Vj are the position and the velocity vectors, respectively. The knowledge of the densities fi^} for all 
n e R 2d and k = 0, 1, 2, defines entirely the state of the system. For a given region f2 the normalization condition 
reads 

00 „ „ 

Mo(*) + y] / drxdvx ... j dr k dv k /4?(n,vi, ... ,r k ,v k ;t) = 1 (2) 
fc=i ^ n J Q 

where /ig (t) is the probability of finding the region void of particles at time t. 

A necessary condition for the occurrence of a pair of particles at the phase space points (r^, vy), (r^v.;) at time 
t > is that Vij —Vi~ Yj , Vij = — Vj belong to the region of the phase space with the characteristic function 

X(r«, vy;*) = dQtij] - a) [l - 0(r y • v«)0 (a - yVd 2 - V* ' ^) 2 ) (3) 




where a denotes the particle diameter, v,j is a unit vector in the direction of the relative velocity, and 9 denotes the 
Heaviside distribution. Indeed, moving backward in time the particles collide during the time interval (0, t) if and 
only if the following three conditions are simultaneously satisfied 

(i) Yij ■ > particles approach each other 

(ii) a > -\/|rjj| 2 — (r^ • v^-) 2 the impact parameter is smaller than a 

(iii) \vij\t > ■ — \j o' 2 — |i"r/| 2 + {y^ ■ ^ij) 2 the time t is long enough 

for the overlapping configuration to occur 

Hence, x( r ij > v ij i ^) = 1) if an d only if no overlapping takes place during the time interval (0,i). 
At time t particles 1, 2, ... k occupy in Q, the one-particles states 

(ri, vi), (r 2 , v 2 ) ... (r fc , v fc ) (4) 

with probability density (0). Using the characteristic function (||) we can construct the probability density for finding 
the same particles in the phase space configuration 

(ri + vidt, vi), (r 2 +v 2 dt,v 2 ) ■■■ (r fc +v k dt,v k ) (5) 
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at time (t + dt), dt > 0. It reads 



Vjjdt, Vjj ; t + dt) 



Mfe (ri,v 1; ... ,r k ,v k ;t) 



In the limit — > + the above expression takes the asymptotic form 



Mfe (n,vi, ... ,r fc ,v fe ;i) 



dt 



dr,. 



x(Tij,Vij-,t)dt 



Using the definition (||) we find 



d_ 

dt 



d 

■fo- ) x(rij,Vij;t) = ■ v y )5(|rij| - a)[l 



where = ry/|ry|. We denote by T v (i,j) the right hand side of (|§|) and rewrite it in the form 

T v (i,j) = a 4 - 1 J da{a • v y )#(-<? • v tf )*(r, 3 - - o£) 



(6) 



(7) 



(8) 



(9) 



Here ? is the unit vector along the line passing through the centers of the spheres at contact. The integration with 
respect to the measure dxr is thus the angular integration over the solid angle in d-dimensional space. The 6 function 
in (JsJ) restricts this angular integral to the hemisphere corresponding to pre-collisional configurations. 
Our aim is to construct the probability density /xj? at time (t + dt) for dt ^ + : 



jti"(ri + vidfjVi, ... ,r k + v k dt,Vk',t + dt) =/4?(r l5 vi, ... r k ,v k ;t) 



(10) 



+ [ §i + J2 v i' -g^rj /4?( r i> v i> -■ r k ,v k ;t)dt 

To this end we have still to add to the term (Q) the probability weights of two events. The first corresponds to the 
presence at time t of (k + 2) particles within in the states 

(ri,vi), (r 2 ,v 2 ), ... ,(r fe ,Vfe),(r + i,v fc+ i),(r fc+2 ,v fc+2 ) (11) 

The state ([ll]) is then transformed into (^J) at time (t + dt) as the result of an annihilating collision between the 
particles (k + 1) and (k + 2), during the time interval (t,t + dt). According to equations (Q) and (||), the rate of 
the occurrence of binary collisions between pairs (i,j) is obtained by applying the operator [— T v (i,j)] defined in (|^) 
to the corresponding distribution. Hence, when dt —> + , the (k + 1, k + 2) annihilation process contributes to the 
density ( |l0| ) the term 

d(k+l) I d(k + 2)T v {k + l,k + 2)^ +2 {l,2 7 ... , k, k + 1, k + 2;t)dt (12) 
n Jn 

where the shorthand notation dj = dvjdvj for j — 1, 2, ... has been used. 

Finally, we have to take into account the effects of the free flow of particles across the boundary dft of the region 
O. Indeed, the fc-particle state can be created or destroyed by an additional particle (k + 1) leaving or entering the 
considered region. Denoting by n the unit vector normal to dil oriented outwards, we get the term 

[dv k+1 [ dS(n-v k+1 )^ +1 (l, ... ,k,k+l;t)dt (13) 
J Jon 

= / d(k + l)v fe+ i • /4? +1 (l, ... ,k,k + l;t)dt 

Here dS is the measure of the surface area, and the equality (Oft follows from Gauss' theorem. 
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The enumerated events combine together to create the complete rate of change of the probability density uP. As 
equivalent events have the same probability measure we can equate (|l0|) with the sum of contributions (f7|) , (p"2|) and 
( |l3| ) obtaining thus the hierarchy equations (k = 1, 2, ...) 



|+X>-^-££nu)K(i, ...,m) = 



(14) 



i< 3 



d{k+l) / d{k+2)T v {k+l,k+2)n k \ 2 (l,2, ...,k+2;t)+ / d(k + l)w k+1 -- 

n Jn Jn or k+1 



Mfc+i(l, - ,k,k+l;t). 



Finally, the evolution equation for p.^ follows from the normalization condition (|^). This completes the derivation of 
the infinite hierarchy of equations satisfied by the probability densities /j,)}. 

From (14) one can derive in a straightforward way the hierarchy satisfied by the reduced distributions 



/fc(l,2, ... , k;t). They are relevant for the evaluation of physical parameters, as //-(1,2, ... , k;t)dl ... dk repre- 
sents the measure of the number of fc-particle phase space configurations, with k particles occupying the one-particle 
states (rj., vi), (r 2 , v 2 ) ... (r^v^) at time t. The distributions f k are related to the probability densities ^ by the 



equation 1 1 1 



A(l,2, ... ,k;t) = f]^±f}l f d(k+l)... f d(k+p)$ +p (l, ... ,k,k+l, ...,k+p;t). 



(15) 



Note that the f k (l, 2, ... , k; t) do not depend on 0. 

In order to derive the evolutio n eq uation for f k one has thus to consider the hierarchy equation ( |l4| ) with k replaced 
by (k + p), and use the relation (15). One finds 



<9r 7 - 



i< 3 



p=0 



P ] - Jn 



d{k+l)... / d(k+ P ) 



(16) 



k+p 

£ 

j=k+i 



k k+p 



k+p k+p 



£ E n<,i)+ E £^) 



j=i i=k+i 



k+l<i <j 



/ife+ p (l, ... ,k+p;t) 



d{k+p+l) / d(k+p+2)T v (k+p+1, k+p+2)p k+p+2 (l, ... , k+p+1, k+p+2;t) 
n Jn 



+ f d(k+p+l) (vfc^x • — — J .» k+p+1; t) 

Jn \ OY k +p+i J 

It is then a question of inspection to see that on the right hand side of (|l^) only the term 

k k+p 



E 

P =i 



(k+p)\ 



j d{k + \)... j d(k+p)J2 E T v (i,j)fi k+p (l, ...,k+p;t) 



j = l i = k+l 



(17) 



f d(k+l)Y,T v (j,k+l)f k+1 (l, ... ,k,k+l;t) 

J 3 = 1 



survives. All the remaining terms exactly cancel out. 

The hierarchy equations satisfied by the reduced distributions f k describing the annihilation dynamics thus read 



A, *-\ A, A.- \ /i A, 

a+E v r^-EE r (M) A(l, ...,k;t) = J d(fc+l)£T"(i,fc+l)/ fe+ i(l, ...,k,k+l;t). (18) 



3=1 



i< j 



3 = 1 
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Consider now equations ( ttq ) supposing that the state of the system is spatially homogeneous. In this case the 
distribution fx does not depend on the particle position. Let us formally take the Grad limit 

(7 — ► 0, n(t) -> oo, ?i(i)CT d_1 = A -1 = const (19) 

where 



n(t) = J dv/i(v;i) 



The fixed mean free path A introduces a relevant length scale, so we pass to dimensionless positions putting 

r j = kKj, j — 1, 2, ... (20) 
With this change of variables the collision operator (||) takes the form 

T v (i, j) = [ n (t)a d ] d - 1 j J da(a ■ v^-S ■ vy)^ - n(t)cx d a] (21) 



We conclude that the term on the left hand side of equation (18) involving the collision operators (|21j) vanishes in 
the Grad limit for d > 1, because the dimensionless parameter n(t)a d tends to zero. Note that this term induces 
dynamical correlations, hindering the propagation of the molecular chaos factorization. On the other hand, using the 
definition of T v (j, k+1) we find that the term on the right hand side equals 



lE [ rfvfc+1 / ' v i(k+i))0(-a- ■ Vj (fcfl) )/ fe+ i(l, ... ,j, ... ,k,r k+ i = Tj-&a,v k +\;t)/n{t) (22) 

3 = 1 3 J 



k 

A 

and its prefactor 1/A remains finite in the same limit. So, formally the hierarchy equations ( |18| ) at a given time t 
reduce pointwise in the Grad limit to the Boltzmann-like hierarchy 

E^'-^T -,*;*) = y d(k + l)Y,T v (j,k + l)f* +1 (l, ...,k,k+l;t), (fc=l,2, ...) (23) 



3=1 3 I " 3=1 



St 

The hierarchy ( [2 3D propagates the factorization of the reduced distributions 

k 



f?(l, ...,k;t) = HfB(j;t) (24) 

3=1 

Hence, if the initial state is factorized, the whole hierarchy ( p3] ) reduces to one non-linear equation 

l +Vl '^) /S(1;t) = /^^(i.^U;*)/ 3 ^*) (25) 



Equation (|2^) is the Boltzmann kinetic equation corresponding to the annihilation dynamics. In the following section, 
we shall see that the formal Grad limit taken here, where n — > oo, is relevant for the description of the annihilation 
dynamics at late times, even if the density n(t) decreases with time. 



B. Scaling analysis of the hierarchy 

The evolution of the annihilation kinetics shares a common feature with the Grad limit: The ratio of particle 
diameter to mean- free-path A = l/(n<7 d_1 ), which is related to the packing fraction, vanishes in both cases. To be 
more specific, we perform a scaling analysis of the exact homogeneous hierarchy equations and look for self-similar 
reduced distributions where the time dependence has been absorbed into the density n(t), with the velocities v and 
positions r renormalized by the typical (root mean squared) velocity v(t) and mean-free-path respectively: We define 
the reduced variables 



v r I f 

— and x = — , with (v) 2 = — — / v 2 fi(v; t) dv. 
v A n(t) J 



(26) 



6 



For the one-body distribution, we therefore introduce the reduced function /i such that 



n(t) 



n{t) 



hiv]t) = W)~ d /l(c) = OT h ^ 



*J(i)° 



v(t) 



(27) 



By definition, the moments of order and 2 of /i are constrained to unity. Requiring that the fc-body distribution 
fk factorizes into Yi^=i m the limit of infinite relative separations between the particles, we consistently obtain 
the scaling form: 



/ fe (ri,vi, ...,T k ,v k ;t) 



/fc(xi,Ci, ...,Xfc,Cfc). 



(28) 



The evolution equations of n(t) and kinetic energy density nv (t) follow from integrating (18) with weights civi and 



v{dvi respectively, for k = 1. We obtain 



dn 

~dt ~ 

d(riv 2 ) 
dt 



= —aw(i) nv 2 , 



where the collision frequency to and kinetic energy dissipation parameter a read 

uj(t) = n(t)v(t) / dc 1 dc 2 d$ (-a ■ c 12 )6(-a ■ c 12 ) f2(c 1 ,c 2 ,a&) 



J d^d^dv (<r • ci 2 ) 6»(-er • c 12 ) c\ f 2 (c 1: c 2 ,aa) 



J c 2 /i(c)dc / dcidc 2 dB (a ■ c 12 ) 6(-a ■ c 12 ) / 2 (ci, c 2 , ctct) 



(29) 
(30) 

(31) 
(32) 



Equation (|3^) is valid for general velocity rescalings; The definition (26) chosen here implies that the term J c 2 /i(c) dc 
in the denominator equals unity. The coefficient a may be seen as the ratio of the kinetic energy dissipated in a typical 
collision normalized by the average kinetic energy, and is time-independent in the scaling regime. It is convenient to 
introduce the internal "clock" C of the dynamics counting the number of collisions, such that dC = udt. With this 
variable, Eqs. (E^) and (E^ol) integrate into 



n{t) = no exp[— C(t)] 



and 



v 2 (t) = v 2 Q exp[-(a-l)C(t)}, 



(33) 



where the time origin with density no and kinetic energy density no«o has been chosen to coincide with C = 0. 
Knowledge of the C-dependence of n and v allows to relate absolute time t to the number of accumulated collisions: 
From Eq. (|3l|), we have 



dC 11 v . , ... 

— = ujq — — = uj exp -C(l + a) 2\, 
dt % Do 



where ujq = u>(t — 0), so that 



The corresponding time evolution is 



C = 



1 + a 



In 1 



1 + a 



■ u t 



n 
no 

v 



1 + a 
2 

1 + a 



co t 



u t 



(34) 



(35) 



(36) 



(37) 



Without knowing the detailed form of the one-particle distribution function /1, it is thus possible to conclude about 
the time decay of n and v, which appear to obey algebraic laws in the long time limit [n(t) cx and v(t) cx t f , 
with £ = 2/(1 + a) and 7 = (a — \)/{a + !)]■ The exponents £ and 7 are consequently simply related to the unknown 
quantity a, for which a perturbative expansion will be put forward in section III before a numerical investigation in 
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section |y|. Moreover, if the initial velocity distribution is of finite support (i.e. vanishes outside a sphere of given 
velocity v*), v fulfills the bound v < v* so that 7 is necessarily positive (or a > 1). In the framework of Boltzmann's 
equation, it will be shown in appendix ^| that the quantity a is necessarily larger than 1. For the specific initial 
condition where all particles have the same kinetic energy at t = 0, v 2 = {v 2 ) is time-independent, and (|3^ ) implies 
that a = 1. From Eq. (|36|), we therefore obtain the time evolution for this situation 

n 1 . 

— = -, 1 38 

n 1 + u> t 

which is exact within the scaling theory. This relation a priori holds in any dimension, except for d = 1 where the 
corresponding initial velocity distribution is the symmetric discrete bimodal function ±c, for which the scaling ansatz 
underlying our approach fails (see the discussion at the end of the present section). 

Inserting the scaling forms (^?J) and ( p8| ) into the first equation of the hierarchy ( |l8| ) imposes the following constraint 
on the decay exponents: £ + 7 = 1. This scaling relation may be simply obtained by elementary dimensional analysis 
@' H' [1 Pi an d ma y be considered as the compatibility condition of the hierarchy with the self-similar scaling 
solutions |12|. It is moreover identically fulfilled by the expressions ( |36| ) and (57). Under the constraint £ + 7 = 1, the 



remaining equations of the hierarchy (k > 1) turn out to be compatible with (28) with the additional information that 
the collision term on the lhs of (|T^ ) decays like i~T~ d ? whereas the remaining terms are associated with a power 1/t. 
Given that 7 + d£ = 1 + (d — 1)£ > 1, this collision term is asymptotically irrelevant except in one dimension where it 
remains of the same order as the dominant ones (1/t). We therefore recover the conclusions obtained by considering 
the formal Grad limit, with distribution functions expected to obey a Boltzmann-like equation. This analysis points 
to the relevance of Boltzmann equation for d > 1, a point which is further corroborated by the numerical results given 
in section |fv|. 

It is interesting to note that both Eq. ( |38| ) and the relation £ + 7 = 1 do not hold in ID for discrete initial velocity 
distributions. In the symmetric situation of a bimodal distribution, the average kinetic energy per particle is conserved 
(so that 7 = 0), whereas the density decays as l/y/t (i.e. £ = 1/2 |Q). The scaling assumption ( p8| ) is consequently 
incorrect in the specific situation of discrete distributions in ID, but valid for continuous distributions ||. To be more 
specific, the scaling form (pq) implies that the collision frequency scales with time like uj cx riv. On the other hand, 
from the analytical solution of the bimodal ±c situation [0 , we obtain oj cx n 2 v with v = c. This discrepancy is the 
signature of dynamical correlations in the latter discrete case. These correlations are responsible for the breakdown 
of feq), and in addition violate molecular chaos. For continuous velocity distributions, again in ID, molecular chaos 



also breaks down while the scaling rt28) is correct. As a consequence, the exponents obtained at the Boltzmann level 



differ from the exact ones (see the discussion in the last paragraph of section III ) , whereas the relation £ + 7 = 1 
holds. 



III. BOLTZMANN KINETIC EQUATION 



This section is devoted to the analysis of the decay dynamics within the molecular chaos framework |b| of the 
homogeneous non-linear Boltzmann equation. No exact solution could be obtained, and our goal is to derive accurate 
approximate predictions for the scaling exponents £ and 7 of the density and root mean squared velocity. 

Before considering the kinetic equation obeyed by the rescaled distribution function, it is instructive to rewrite the 
original homogeneous Boltzmann equation (|2q) in the form 



dfi(y;t) 

at 



-z/(v;i)/i(v;i) with f(vi;t) 



a d 1 J da{a ■ Vi 2 )0(<7 • v 12 ) / dv 2 |vi 2 |/i(v 2 ; t) 



(39) 



where in the last equation, the term in brackets may be computed explicitly as function of dimension (it is understood 
that the unit vector v 12 denotes an arbitrary direction). For our purpose, it is sufficient to notice that at all times, 
the collision frequency z/(v; t) of the population having velocity v remains finite in the limit v — > 0, provided the 
first moment of fx exists. In this situation, Eq. ( |39"|) implies that /i(v;t)//i(v;0) admits a finite limit for v — > 0, 
or equivalently, it may be stated that if the initial velocity distribution behaves like i> M near the velocity origin, this 
feature is preserved by the Boltzmann dynamics. Previous works have shown accordingly that the scaling exponents 
£ and 7 depend on the exponent ^ [|, [l(J . 

Making use of relations (^) and (po|), insertion of the scaling form ( p7| ) into the Boltzmann equation leads to 



1 



1 - a 



d 

d + ci — 

dc\ 



/i(ci) = fxict) 



dc 2 -p^r /i(c 2 ), 

C12 



(40) 
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where we have assumed an isotropic velocity distribution [/i(c) = /i(c)] and introduced the average ((■■•)) = 
/(. . .)/i(ci) f\{c-i)dc\dc2- (C12) is therefore the rescaled collision frequency. Once fi has been chosen, Eq. (|4(i| ) 
admits a solution for a unique value of a. We show in appendix |A| that the inequality a > 1 necessarily holds. 

Irrespective of a, the large velocity behaviour of /1 may be obtained following similar lines as in ||. pL [fi], |l5): it 
is possible to integrate formally (Eoh and cast /1 into 



/i(c) 
/i(c') 



(7) 



exp 



1 - a (ci2 



?(c") 



dc" 



In this equation, j/ is itself a functional of /1 : 



C12 /i(c2) dc 2 , 



such that v(c)jc goes to a finite limit for c — * 00. We therefore obtain the large velocity tail 

2 c 



~ 2 + d(l- Q ) 

/i(c)occ 



exp 



1 c 12 



for 



(41) 



(42) 



(43) 



In one dimension, we recover the results of references Q] and (^] . In (J, |E| , an approximation was derived for a [or 
equivalently £ = 2/(1 + a)] by assuming that the large velocity behaviour of /1 could hold for the whole velocity 
spectrum. In this picture, the power of c appearing on the rhs of Eq. ( [i"3| ) is equated to the exponent /i characteristic 
of the small velocity behaviour (imposed by the initial distribution chosen, see above), with the result 



a 



2d + 2fi 
2d + 2ji+ 1' 



(44) 



This prediction encodes the correct dependence on [i and dimension (£ increases when /1 or d increase) , and turns out 
to have an accuracy of order 10% when compared to the numerical results [fl0|| . In the limit of large dimension, we 
obtain from (Q) £ ~ 1 — (2d) -1 , whereas Krapivsky and Sire have shown that £ ~ 1 — dr x (\ — l/\/2), also in the 
framework of the Boltzmann equation. The remainder of this section is devoted to the derivation of a more precise 
value for a, which furthermore coincides with the exact 1/d correction for d — * 00. 
Invoking the identity 



do-" ( d + c- )/i(c) 



-p {c- 



the energy dissipation parameter a may be given the set of equivalent expressions: 

2 
P 



1 



(Cl2 C\ ) 



1 



Vp > 0. 



A particularly useful relation between a and moments of fx follows from considering the limit c\ — > of 

2 A ( Cl ) 



a 



1 



1 



C12 



(45) 

(46) 
we get 
(47) 



The (infinite) family of relations ( |47| ) and (^6| ) is equivalent to the original integro-differential equation (fio]), and 
well suited to a perturbative analysis. To this end, a systematic approximation of the isotropic function f\ can be 
found by expanding it in a set of Sonine polynomials |]l6| : 



h(c) = M(c) 



1 + fln^^C 2 ) 



(48) 



These polynomials are orthogonal with respect to the Gaussian weight 



M(c) = ( - 



d/2 



,-dc 2 /2 



(49) 
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and the first few read 



S (x) = 1 (50) 
d 
2 



S ± (x) = ^{-x + 1) (51) 



= (52) 
o 4 o 

The coefficients a n follow from the orthogonality relation J S n {c 2 )S m {c 2 )M.(c)dc oc S nm : In particular, 

ai-^(^i( C 2 )> = ^(l-( C 2 »=0 (53) 



from the definition of rescaled velocities fl2q). The first non-Gaussian correction is thus embodied in ci2, that is 
proportional to the fourth cumulant (kurtosis) of the velocity distribution: 

oa = y [(4>-3( C 2 ) 2 ] =^2^)-!, (54) 

with Cj a Cartesian component of c. Upon truncating Eq. ( ff8| ) at a finite order n, we obtain a regular velocity 
distribution near c = 0. We consequently restrict our analysis to the case fi = (the dependence on /i has been 
considered in pj). 

It is also noteworthy that any truncation of ( [48| ) at arbitrary order n leads to a Gaussian high energy behaviour, 
incompatible with the resu lt (f43| ) corresponding to an overpopulated tail with respect to the Maxwellian. However, 
it will be shown in section |lV[ that the difference between the truncated expansion ( |48| ) and the numerical velocity 
distribution becomes manifest far in the tail, where the distribution has reached very low probabilities. Consequently, 
when the moment involved in ( [46| ) and (E^ are evaluated from the truncation of (f48|), the accuracy of the result 
is expected to be better for low orders p in ([f6|). Hence the privileged role played by (|47|), which is of lower order 
than any of the identities (f46|). In practice, upon truncating (48) at order n, the n unknowns a, a 2 , ■ ■ ■ , a n are 
computed evaluating the corresponding moments appearing in ( 47 ) and in n — 1 of the relations (J46|) , among which 
it is convenient to retain the n— 1 even values of p (p — excluded). Truncation of ( [48| ) at order n — 2 yields precise 
predictions for a and /i, and already at Gaussian order, a turns out to be very close to its numerical counterpart: 
Setting n = (or equivalently n = 1 since a± = 0) in (|48|), we obtain immediately the zeroth order approximation 




a = a = 1 + - j I — ! (.V-,, 

which corresponds to 

2 2d 

^ e ° = TT^ = 2(d + i)-V2- (56) 

For large dimensions, this estimation goes to unity, with the exact l/d correction computed in H 

(57) 




This behaviour turns out to be "universal" , in the sense that the \i dependence does not appear at this order [jl0[ . 

We shall also be interested in the non-Gaussian features of the velocity distribution, that we quantify by the fourth 
cumulant 02- The (cumbersome) calculations at second order in Sonine expansion are detailed Appendix [b|. We 
obtain 

q 2 =8 < 2 ^- 3 ) (58 ) 
4d 2 + d(6-72) V ; 

V2 

a 2 = a + —— a 2 . (59) 
loot 

The corresponding density exponent follows from £ 2 = 1/(1 + a 2 ) as before. For d — > 00, et 2 ~ 2(3 — 2^/2)d~ 1 , 
irrespective of fi []10[| which reinforces the "universal" nature of large d. The correction to £ carried by a 2 behaves 
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as 1/d 2 in this limit, and does not affect the X/d terms which still coincide with the exact behaviour ( pT\l. Both 
predictions (55) and (|59| ) are such that a > 1, which is required to obtain a normalizable distribution in (]4l]) and 

The second order expansion considered here may be improved by consideration of higher order Sonine terms and 
inclusion of non- linear terms in et2. In the related context of inelastic hard spheres, the limitation of working at linear 
order in ai with neglect of Sonine terms of order n > 3 may be found in |l7|. Alternatively, keeping non linear terms 
in ci2 and neglecting again Sonine terms of order n > 3 leads to multiple solutions. A stability analysis is then required 
to determine which one is stable, as discussed in fl8|| . 

Here, the value obtained for a 2 is quite small (see below). Our approximate expressions are accurate when compared 
to the full numerical solution of Boltzmann's equation, so that we did not calculate any higher order coefficients, nor 
consider non linear terms in 02- The existing body of literature reports, within Boltzmann framework, numerical 
exponents in ID that are in excellent agreement with our predictions, already at zeroth order. For the case /1 = 0, 
expressions (|56|) and ( |59| ) give £ — 0.773 and £2 — 0.769 whereas the numerical result obtained in || is £ ~ 0.769. 
These exponents are close to their counterparts extracted numerically from the exact dynamics (0.785 ±0.005 in ||, 
and more recently 0.804 [^9|). The difference between the exact exponents and those obtained assuming molecular 
chaos is consistent with the conclusion of section |0]: In 1 dimension, the factorization of the two-body distribution 
fi underlying Boltzmann ansatz is not an exact property of the distributions obeying hierarchy (p^|). On the other 
hand, for d > 1, the molecular chaos exponents are expected to become exact. This property is illustrated in the next 
section. 



IV. SIMULATION RESULTS 



The numerical results presented in this section correspond to the situation /i = 0, unless stated. We refer to |T(| 
for the case of diverging (fi < 0) or depleted (/z > 0) velocity distributions near v = 0. 



A. The numerical methods 



We follow two complementary numerical routes. First, we solve the time-dependent homogeneous Boltzmann 
equation by means of the Direct Simulation Monte Carlo method (DSMC), originally developed to study ordinary 
gases ptj. This scheme, where a suitable Markov chain is constructed, has been extended to deal with inelastic 
collisions [^l], and is easily modified to describe the situation under study here, which does not conserve the total 
number of particles. Restricting to a spatially homogeneous system, the algorithm is especially easy to implement, 
and may be summarized as follows: among Nq initial particles having a given velocity distribution, a pair (i,j) is 
chosen at random, and removed from the system with a probability proportional to |vy|. The (suitably renormalized) 
time variable is then incremented by the amount (N \vij\) , where N is the number of particles remaining in the 
system, before another pair is drawn etc... This scheme provides the numerical exact solution of Eq. (|39|), and 
allows to test the validity of the approximations put forward in section III [essentially: Truncation at second order in 
expansion ( f48| ) supplemented with calculations performed at linear order in the fourth cumulant 02, see Appendix |b| . 
A precise analysis of the late time dynamics (and especially the computation of velocity distributions) suffers from 
the concomitant low number of particles left, and the statistical accuracy is improved by averaging over independent 
realizations. 

The second numerical method (Molecular Dynamics p3| ) consists of integrating the exact equations of motion for 
an assembly of spheres confined in a (hyper)-cubic box with periodic boundary conditions. This route assesses the 
validity of the approach relying on the homogeneous Boltzmann equation, but does not offer the same accuracy as 
DSMC, nor the possibility to follow the evolution over comparable times. In particular, once the mean-free-path A, 
which grows rapidly as Ifi, becomes of the order of the box size L, the subsequent evolution suffers from finite size 
effects and should be discarded: When A > L, the algorithm used is unable to find collision events for those particles 
which make more than one free flight round on the torus topologically equivalent to the simulation box, which causes 
a spurious slowing down of the dynamics. It is then tempting to reduce A by increasing particles diameter a, but 
then finite (and necessarily transient) density effects -incompatible with the scaling assumption (28)- may also arise 
if the packing fraction 0, proportional to na d , is not low enough. Simulating explicitly the limit of point particles, 
the DSMC scheme considered here is free of this defect. The initial number of particles considered in MD needs to 
be large to allow the system to enter the scaling regime before finite-size effects become dominant: we considered 
systems with N = 5.10 4 to N = 5.10 5 spheres initially (compared to N — 10 6 to 10 s in DSMC). 
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B. Dynamic scaling behaviour 

The results of two-dimensional DSMC and MD simulations are shown in Figure [l], where the quantity on the x 
axis is expected to scale like real time t from the scaling relation £ + 7 = 1 . This log- log plot is a direct probe of the 
exponent £, from the slope measured. Both MD and DSMC methods give compatible results, with the possibility to 
follow the dynamics over a longer time interval in DSMC. The departure observed for rio/n ~ 200 corresponds to the 
slowing down of MD dynamics resulting from finite-size effects (see Figure || below). The theoretical predictions at 
zeroth and second order are very close (£0 = 0.872 and £2 = 0.870), and in excellent agreement with the simulation 
results over several decades. On a similar graph as Fig. frl the kinetic "temperature" v 2 exhibits a lower law behaviour 
(not shown) with an exponent —2j in good agreement with the theoretical prediction (7 ~ 0.13 in 2D). Moreover, 
the exponents obtained analytically and numerically are compatible to those reported in the literature (£ ~ 0.89 in 
H and £ ~ 0.87(5) using a multi-particle lattice gas method pi|). 

The time evolution of inverse density and inverse typical velocity square is shown in Figure ^, where considering 
no/n — 1 and v 2 /v 2 — 1 instead of ixq/ti and Vq/v 2 allows to probe the short time behaviour. Unless stated, the 



initial velocity distribution is an isotropic Gaussian. From Eqs. (29) and (J30|), n and v 2 evolve linearly with t for 



uj t <C 1 [see also ([36|) and (37)]; the same holds for inverse density and inverse typical velocity squared, which is 
indeed observed in Fig. ||. MD and DSMC result super-impose, except at late times where MD dynamics suffers from 
the slowing down discussed previously. For both numerical methods, the scaling relation £ + 7 = 1 is well obeyed, in 
principle at late times only, in the scaling regime. Special combinations of n and v can however be constructed with 
the requirement to match the short time evolution with the scaling behaviour. One of these quantities is displayed in 
Fig. ||, with a resulting scaling regime extending over more than 10 decades in time. In Fig. ^, we not only test the 
validity of the theoretical scaling exponents, but also the full time dependence as predicted by Eqs. (|36|) and (|37|). In 
order to improve the agreement between theory and simulation (which holds over more than 6 orders of magnitude in 
time), the system has been left time to enter the scaling regime: The time origin t = has been chosen when 80% of 
the particles originally present have disappeared. The corresponding reference configuration (with subscripts 0) thus 
differs from the ones considered previously. 



C. Velocity distribution in the scaling regime 

In order to understand the reasons for the good agreement between our theoretical predictions and the simulations, 
we now consider the velocity distribution, restricting to Monte Carlo results (leading to similar ^conclusions, MD 
is much more demanding on CPU time and does not allow to investigate detailed features of f± with the same 
accuracy). After a transient where the probability distribution function fi evolves with time, a well defined scaling 
regime is reached with a time independent f± (c) shown in Fig. together with the Sonine prediction pushed at second 
order. The agreement is remarkable, and it is also observed that the Gaussian approximation is already close to the 
asymptotic rescaled velocity distribution, which is quite surprising in a kinetic process extremely far from equilibrium, 
with furthermore no conservation laws. Given that our perturbative analytical work relies on the calculation on low 
order moments of /1, this explains the accuracy of the zeroth order predictions o.q and £o- From Eq. J43|), we expect 
the differences between the Sonine expansion and the numerical distribution to become visible in the high energy tail, 
which is confirmed by Fig. |(| As predicted, f\ is overpopulated with respect to the Gaussian, and displays a high 
velocity tail of the form (|4^) (see the inset). 



D. Evolution toward the asymptotic solutions 

Before the scaling regime is attained, /1 is time-dependent, as shown in Fig. |7], where the distributions at different 
times have been renormalized by M. to emphasize the building-up of non-Gaussianities. The evolution towards the 
scaling solution 1 + 02152 can be observed. With respect to the Gaussian, f\ is at all times overpopulated both at large 
and small velocities (which may be related to the positive sign of 02 for the latter case) ; normalization is ensured by 
an under-population at intermediate velocities. Figures |5|, |^ and [?] show that the indirect measure of 02 through the 
non-Gaussian character of fi/M agrees with the theoretical prediction, but it is also possible to compute directly 02 
in the simulations through its definition as a fourth cumulant [Eq. (|54j)]. It turns that both methods are numerically 
fully compatible. Moreover, the Sonine expansion (|48|) truncated at n = 2 holds at any time, even in the transient 
regime, with the time dependent fourth cumulant 02 measured from ([34]) (see Fig. ||). This result is not a priori 
expected and points to the relevance of the expansion jiq) . We did not try to solve analytically the homogeneous 
time-dependent Boltzmann equation within the same framework as in the scaling regime, so that we do not have any 
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predi ction for the (transient) time-dependence of 0,2 and a. However, as shown in the inset of Fig. the relation 
(plOj) (which reads in 2D a = 5/4 + 7a2/16) remarkably holds for all time. Here, the energy dissipation parameter 
a has been computed through the ratio (C12 c\) / {(c\2) {c\)) = (cf ) co n/(cf ) , where (cf) co n, the mean energy dissipated 
in a collision is computed in the simulations and normalized by the time-dependent mean kinetic energy per particle 
(cf). 

E. A final remark: Identification of "isobestic" points 

For /i = 0, Figure |?j indicates that during the transient evolution toward scaling, the distributions of reduced 
velocities have fixed points (for a given initial velocity distribution, the curves corresponding to /1 at different times 
all pass through common points, that we called isobestic points). This feature has been observed for all initial 
distributions investigated and appears to be a systematic property of the dynamics, which still holds for non vanishing 
values of \x (see Fig. ^|). We have no analytical explanation for this observation. 

V. CONCLUSION 

An analytical derivation of the equations governing the dynamics of an infinite system of spherical particles in a 
c?-dimensional space, moving freely between collisions and annihilating in pairs when meeting, has been obtained. 
The hierarchy equations obeyed by the reduced distributions /fe(l, 2, A;; t) have been derived. In the Grad limit, 
this hierarchy formally reduces to a Boltzmann-like hierarchy for dimensions d > 1. If these reduced distributions 
fu factorize at the initial time, this factorization remains valid for all times and the whole hierarchy reduces to one 
non-linear equation. 

In the long time limit, the ratio of the particle radius to the mean-free path vanishes. A scaling analysis of the 
exact homogeneous hierarchy equations has been performed. Self-similar reduced distributions in which the time 
dependence has been absorbed into the density n(t) and the root mean-square velocity v(t) were introduced. As a 
result, the exponents £ and 7 describing the decay with time of n(t) and v(t) depend only upon one single parameter 
a, related to the dissipation of energy. Moreover, it turns out that in dimension higher than 1, the terms responsible 
for the violation of molecular chaos are asymptotically irrelevant. Therefore, we recover the conclusions reached in the 
formal Grad limit and thus, the Boltzmann equation becomes exact in the long time limit in dimensions higher than 
1. The above arguments give a first principle justification for the use of the Boltzmann equation approach for d > 1, 
as well as its limitations, in ballistic annihilation problems, an issue that has been overlooked so far. However, as 
discussed above, this scaling analysis is incorrect for a 1 dimensional system with discrete initial velocity distribution, 
a situation for which the scaling assumption underlying our approach fails. 

The Boltzmann equation has been solved within an expansion in Sonine polynomials S n . Truncation to order n = 2, 
provides the first non-Gaussian corrections to the scaled velocity distribution, and leads to analytical predictions for 
the exponents £ and 7 as a function of d. For large dimension d, these predictions coincide with the exact l/d 
correction to the naive mean field values (£ = 1 and 7 = 0), calculated in Q. The above analytical predictions are 
in remarkable agreement with the results of extensive numerical simulations we have performed for two dimensional 
systems (implementing the complementary Monte Carlo and Molecular Dynamics techniques). In ID for regular 
continuous velocity distributions, it is noteworthy that they are in excellent agreement with the numerical solution of 
the Boltzmann equation, and quite close to the exact values obtained with Molecular Dynamics (4% difference). This 
last point was unexpected since molecular chaos breaks in ID. Finally, the time dependence of the reduced velocity 
distribution function shows an unexplained and remarkable feature, with the existence of fixed ("isobestic") points, 
irrespective of initial conditions. 
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APPENDIX A 

Within Boltzmann's kinetic equation, we show in this appendix that a > 1. To this end, Eq. ( pp[ ) is rewritten in 
the form 

Uci) + (^) ^ ( C l/l( C l)) = /l( C l) / dc 2 7l(C2). ( A1 ) 

and integrated with weight ijj(ci)dci, where we choose ip(c) — c 2 ^ x ^ a \ Assuming that the moment J ipfi exists, we 
obtain after an integration by parts (neglecting surface terms) 

0= [dcidc^{cx)-^-}i{ci)Jx{c2). (A2) 

J (Cl2> 



The right hand side of Eq. (A2) is a strictly positive quantity [except when /i(c) = 5(c)], which leads to a contradiction. 
We therefore conclude that the quantity J ipf\ does not exist. Remembering that -0(c) /i(c) ~ c ^+ 2 /( 1 ^ Q ) near the 
velocity origin, the divergence of J ipfi implies 

o 

d-l+n<-l. (A3) 



1 - a 



Supplementing this condition with the normalization constraint /i > — d, we obtain a > 1. The physical meaning of 
this condition is that the typical kinetic energy dissipated per particle in a collision is larger than the average kinetic 
energy of the system at the same time. The temperature v 2 is therefore a decreasing function of time. 

APPENDIX B: SECOND ORDER TRUNCATED SONINE EXPANSION 

In this Appendix, we calculate the dominant deviation of fi from the Gaussian shape, and the associated energy 
dissipation parameter a by setting /i(c) = M(c) {l + a2<S2(c 2 )}. The moments appearing in ( [47| ) and ( f4(f ) with p = 2 
are then computed as a function of a 2 . Upon writing 

(ci) \ (ci 2 cf) 



a = 1 + — =- 1- f^- = / "J" , (Bl) 

the last equality provides an equation for 02 which is solved so that a is finally explicitly known as a function of input 
parameters /1 and dimension d. As in the main text, the angular brackets denote averages with weight /i(ci)/i(c2) 

(■■■)fi(ci)fi(c 2 )d Cl dc 2 (B2) 
(...)M(ci)A4(c 2 ){l + a 2 [S 2 (cl) + S 2 (c 2 2 )]} d Cl dc 2 + O (af) . (B3) 

In the following, non-linear terms of order a 2 will be neglected. In the spirit of ^5|, it is convenient to introduce 
center-of-mass and relative velocities Ci = C + c\ 2 /2\ c 2 = C Ci 2 /2 and to define the generic moments 

M np = (c n 12 CP) (B4) 



/ da 2 dCc 12 C p I — I e 12/ 1 1 + 02 y(ci+c 2 ) (c x + c 2 J + > 



(B5) 

From cf + c 2 — 2C 4 + 2(C-Ci 2 ) 2 + c\ 2 /S + C 2 c\ 2 , the term (C-Ci2) 2 appearing under the integral sign in ( |B5| ) becomes 
C 2 c\ 2 /d for symmetry reasons, and making use of c\ + c\ = 2C 2 + c\ 2 /2, the variables c 12 and C decouple in QB5| ). 
The resulting integrals yield 

^ = = 1 + $d {d{n2 + p2) 2d{n +p)+ 2np(d + 2)} + {al) ■ (B6) 

In this equation, the subscript refers to averages with Gaussian measure (formally a 2 = 0): 

r ( d+n ) r ( d+p } 

(c? 2 C p >o = (Vd)— p 2" 2 r2(d/2 V } 2 7 , (B7) 
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where T is the Euler function. The moments (C12) and (ci2 cf) — (ci2(c 2 + c 2 ,)) /2 = (cf 2 )/4 + (C12 C 2 ) appearing in 



(Bl) are then known 



(C12) = 

(Cl2) 
(fl2 C 2 ) 

(Cl2> 
(C12C 2 ) 

(ei2)(4) 



16 08 



0(a 2 ) 



= 1 



1 

2d 
1 

2d 



ft 2 

IT 



0{al) 



(B8) 
(B9) 
(BIO) 



For an elastic hard sphere fluid at equilibrium (with thus 02 = 0), this last quantity equals 1 + l/(2d) and represents 
the ratio of the mean kinetic energy of colliding particles (averaged over successive collision events) to the mean kinetic 
energy of the population. As expected, this ratio exceeds 1, since typical colliding partners are "hotter" than the 
mean background. This quantity is easily measured in Molecular Dynamics or Monte Carlo simulations (see e.g. Fig. 
|J). We also note that the ratio (B8) has been computed in |Q at the same level of approximation in the context of 
rapid granular flows, with the same result, van Noije and Ernst also reported non Gaussian corrections to the cooling 
rate 7 of an inelastic hard sphere fluid 113] 



7o 



1 + —a 2 + O (aj) Vd, 



(Bll) 



where 70 denotes t he co oling rate evaluated assuming Maxwellian velocity distributions. In terms of the mom ent s 



M np introd uced in (B4), it can be shown that 7/70 = M^o/M^ and it is then easily checked that expression (B6) 



reduces to (Bll) for n — 3 and p = 0, for arbitrary dimensionality. 

The remaining unknown quantity is (ci) which may be calculated following similar lines as above: 



1 



'1/0 



'1/0 



0-2 1 

—n(n 



^/2 



2) + O (a 2 ) 



r(|) 



from which we extract (ci) = (c). As expected, the 0,2 correction in (|B13 |) vanishes for n = 
respectively from the normalization constraint and the definition (|26|) of c implying (c 2 ) 
obtain @ and (§§) from Eqs. (pi). 



(B12) 
(B13) 

and n = 2, which follows 
= 1. Gathering results, we 
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FIG. 1: Inverse density no/n versus noVo/(riv), for d — 2. The dotted line has slope £2 = 0.87 [see equation (Q6|)]. Initial 
number of particles: 5.10 6 (with a further average over 10 3 independent replicas) for DSMC and 2.10 5 for MD. In both cases, 
the initial velocity distribution is Gaussian (fi = 0), and the initial configuration used for MD is that of an equilibrium hard 
disk fluid with packing fraction = 0.1 (chosen low enough to avoid finite packing effects). 



10 




FIG. 2: Plots of (no/n) — 1 [upper curve corresponding to DSMC, compared to its MD counterpart (crosses)] and (vo/v) 2 — 1 
(lower dashed curve for DSMC, circles for MD), as a function of time. The dotted line at short times has slope 1. 

15.0 i . 1 . 1 1 




-5.0 0.0 5.0 

log 10 ((Q t) 

FIG. 3: Plot of log 10 [no/n — 1] + log 10 \noVo/(riv 2 ) — ll on a logarithmic time scale (d = 2). DSMC results are represented by 
the continuous curve, while the crosses correspond to MD. The dashed line has slope 2. At late times, the quantity displayed 
is expected to behave as 2 log 10 (^ot)_from the scaling relation £ + 7=1. At short times, the same behaviour is observed, for a 
different reason [see Eqs. (M) and (B7f)] ■ The ultimate MD slowing down is again visible. 
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ay 

FIG. 4: Time dependence of n (lower curve) and v 2 (upper curve) obtained in Monte Carlo, compared to the dashed curves 
corresponding to the theoretical predictions (^) and (p7|) , where the energy dissipation coefficient a is calculated at second 
order in Sonine expansion [012 — 1.297 from Eq. (|59[)]. 




1 ■ 1 ■ 1 ■ 1 

0.5 1 1.5 

C , 

FIG. 5: Probability distribution function /i(cj) of a given Cartesian component Ci of the rescaled two-dimensional velocity c. 
The time independent distribution obtained in DSMC simulations at late times is compared to the Gaussian M and the Sonine 
expansion truncated at n = 2, with the fourth cumulant given by Eq. (^) (0,2 — 0.109 for d = 2 and y, = 0). All distributions 
have variance 1/2 so that (c 2 ) = 1. The results have been obtained by averaging over 50 replicas of a system with initially 
TV = 40. 10 6 particles. 
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C. 

i 

FIG. 6: Same as Figure |B|, on a linear-log scale to probe the tails of the distributions. In the inset showing c~ 4 ' 7 fi(ci) as a 
function of Cj, the dashed curve corresponds to the Gaussian [c~ 4 ' 7 M(a)] while the straight line is a guide for the eye evidencing 
a pure exponential behaviour. For d — 2 and fi — 0, the exponent d + 2/(1 — a) appearing in Eq. ((43) is close to -4.7 and has 
been used to rescale the quantity plotted on the j/-axis in the inset. 




0.85 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

-3.5 -2.5 -1.5 -0.5 0.5 1.5 2.5 3.5 

C. 

i 

FIG. 7: Plots of f(ci)/M{ci) versus d, at different times corresponding to the indicated densities. The initial distribution is 
Gaussian (thus corresponding to the flat curve), and the thick curve is Sonine solution 1 + a252 with ai given by (p8|). 
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FIG. 8: Plot of f(d)/A4(ci) as a function of Ci, at the particular time t 1 ^ 2 where the density is exactly half the initial one, 
with the same initial distribution as in Fig. [?] [the circles (DSMC measure) thus show the same distribution as the circles of 
Fig. fjj]. The thick curve shows 1 + £12(^/2) £2 with a^itx/^) measured from its definition ([Hi]). The dashed curve is the Sonine 
prediction in the scaling regime (i.e. the thick curve of Fig. Q). Inset: a as a function of 0,2 (see main text). The DSMC 
measure is compared to the prediction ( |B10| ) shown by the straight line ending at the point -indicated by a cross- of coordinates 
(0.109,1.297) as predicted by Eqs. (^) and (|59|). The square located at (0.0207,1.2595) corresponds to the numerical measure 
of a and a 2 made at time ti/2 for which the velocity distribution is displayed in the main graph. 





FIG. 9: Plots of fi(ci) as a function of Ci at different times. The left graph corresponds to an initial velocity distribution with 
fi = —3/2 while /i = 3 for the right graph. On both graphs, the initial distribution is shown by the dashed curve, whereas the 
thick curves display the asymptotic distributions approached in the scaling regime. 



